#!/usr/bin/env Rscript

# 02_clustering_dim_reduction.R
# Seurat聚类和降维分析
# 输入: seurat_obj_qc.rds
# 输出: seurat_obj_clustered.rds, 多种可视化图形

library(Seurat)
library(ggplot2)
library(cowplot)
library(dplyr)

# 设置随机种子保证结果可重复
set.seed(123)

cat("=== Seurat聚类和降维分析开始 ===\n")

# 读取QC后的数据
cat("读取QC后的Seurat对象...\n")
seurat_obj <- readRDS("output/seurat_obj_qc.rds")

# 数据标准化
cat("执行数据标准化...\n")
seurat_obj <- NormalizeData(seurat_obj, 
                           normalization.method = "LogNormalize", 
                           scale.factor = 10000)

# 识别高变基因
cat("识别高变基因...\n")
seurat_obj <- FindVariableFeatures(seurat_obj, 
                                  selection.method = "vst", 
                                  nfeatures = 2000)

# 提取前10个高变基因
top10 <- head(VariableFeatures(seurat_obj), 10)

# 绘制高变基因图
variable_plot <- VariableFeaturePlot(seurat_obj)
variable_plot <- LabelPoints(plot = variable_plot, 
                            points = top10, repel = TRUE)

# 保存高变基因图
ggsave("output/variable_features.png", variable_plot, 
       width = 10, height = 6, dpi = 300)

# 数据缩放
cat("执行数据缩放...\n")
all_genes <- rownames(seurat_obj)
seurat_obj <- ScaleData(seurat_obj, features = all_genes)

# 线性降维（PCA）
cat("执行PCA降维...\n")
seurat_obj <- RunPCA(seurat_obj, features = VariableFeatures(object = seurat_obj))

# 可视化PCA结果
cat("生成PCA可视化...\n")

# PCA特征图
pca_features <- VizDimLoadings(seurat_obj, dims = 1:2, reduction = "pca")
ggsave("output/pca_loadings.png", pca_features, width = 10, height = 6, dpi = 300)

# PCA散点图
pca_plot <- DimPlot(seurat_obj, reduction = "pca")
ggsave("output/pca_plot.png", pca_plot, width = 8, height = 6, dpi = 300)

# PCA热图
pca_heatmap <- DimHeatmap(seurat_obj, dims = 1:6, cells = 500, balanced = TRUE)
ggsave("output/pca_heatmap.png", pca_heatmap, width = 10, height = 12, dpi = 300)

# 确定数据集的维度
cat("确定数据维度...\n")
pca_elbow <- ElbowPlot(seurat_obj)
ggsave("output/pca_elbow.png", pca_elbow, width = 8, height = 6, dpi = 300)

# 细胞聚类
cat("执行细胞聚类...\n")
seurat_obj <- FindNeighbors(seurat_obj, dims = 1:10)
seurat_obj <- FindClusters(seurat_obj, resolution = 0.5)

# 非线性降维（UMAP/t-SNE）
cat("执行非线性降维...\n")

# UMAP
seurat_obj <- RunUMAP(seurat_obj, dims = 1:10)

# t-SNE
seurat_obj <- RunTSNE(seurat_obj, dims = 1:10)

# 可视化聚类结果
cat("生成聚类可视化...\n")

# UMAP聚类图
umap_clusters <- DimPlot(seurat_obj, reduction = "umap", label = TRUE)
ggsave("output/umap_clusters.png", umap_clusters, width = 8, height = 6, dpi = 300)

# t-SNE聚类图
tsne_clusters <- DimPlot(seurat_obj, reduction = "tsne", label = TRUE)
ggsave("output/tsne_clusters.png", tsne_clusters, width = 8, height = 6, dpi = 300)

# 保存聚类后的对象
cat("保存聚类后的Seurat对象...\n")
saveRDS(seurat_obj, "output/seurat_obj_clustered.rds")

# 生成聚类统计信息
cluster_stats <- table(seurat_obj@meta.data$seurat_clusters)
cat("聚类统计:\n")
print(cluster_stats)

# 保存聚类统计
write.csv(data.frame(Cluster = names(cluster_stats), 
                    CellCount = as.numeric(cluster_stats)),
          "output/cluster_statistics.csv", row.names = FALSE)

cat("=== 聚类和降维分析完成 ===\n")
cat("输出文件:\n")
cat("- output/seurat_obj_clustered.rds: 聚类后的Seurat对象\n")
cat("- output/variable_features.png: 高变基因图\n")
cat("- output/pca_*.png: PCA相关可视化\n")
cat("- output/umap_clusters.png: UMAP聚类图\n")
cat("- output/tsne_clusters.png: t-SNE聚类图\n")
cat("- output/cluster_statistics.csv: 聚类统计\n")